Grazing effects are multi-trophic through mediating plant composition: a meta-analysis



Alessandro Filazzola, Batbaatar Amgaa, Charlotte Brown, Issac Heida, Jessica Grenke, Margarete Dettlaff, Tan Bao, & JC Cahill

library(tidyverse)
library(PRISMAstatement)

Purpose

Conduct a meta-analysis of the literature testing the indirect effects of grazing on animal taxa’s through the direct effects on the plant community.

Objectives

Timeline

date task
Nov 9 Establish search terms to be used in the meta-analysis
Nov 12 Compile list of journal artcles and sub-divide for each researcher
Nov 14 Begin reviewing papers and extracting data
Jan 28 Complete data extraction from papers
Feb 11 Complete preliminary analysis and set structure for MS
Feb 25 Settle on analyses to be used and begin writing manuscript
March 11 Complete first draft of MS and pass to co-authors
March 25 Comments passed back on draft
April 2 Complete revisions and submit to journal

Literature Review - 1. Search

Search Terms

A systematic literature search was conducted using Web of Science for all emperical research articles. The review will include all studies globally. The intended purpose of this search is to capture all articles that have documented grazing either along a gradient (e.g. different frequencies or intensity) or presence/absence (e.g. excluded, ungrazed, or retired ranch lands). We condcuted two separate searches to capture studies that tested gradients and studies that compared grazing to ungrazed treatments. Duplicate articles between the searches were removed. We also intentionally excluded terms that resulted in articles not relevant to the purpose of this study including: review, synthesis, policy, social, carbon, and fish. The search terms that used were:

Search A graz* OR livestock AND exclosure* OR exclusion OR exclude* OR ungrazed OR retire* OR fallow*

Search B grazing intensity OR grazing gradient OR stocking rate OR rotation* grazing

Literature Review - 2. Sort

This steps includes a. checking for duplicating, b. reviewing each instance for relevancy, c. consistently identifying and documenting exclusion criteria. Outcomes include a list of publications to be used for synthesis, a library of pdfs, and a PRISMA report to ensure the worflow is transparent and reproducible. Papers were excluded with the following characteristics:

  • Not emperical study (e.g. review, book chapter)
  • Irrelevant categories (e.g. political science, law, sports tourism, art)

Reasons for exclusion

evidence <- read.csv("data//ReviewerSearch//evidence.new.csv")
### Identify studies that were excluded
excludes <- evidence %>% group_by(reason.simplified) %>% count(reason.simplified) %>% filter(reason.simplified!="")

## Generate plot
excludes$reason.simplified <- factor(excludes$reason.simplified, levels = excludes$reason.simplified[order(excludes$n)]) ## sort by frequency

ggplot(excludes %>% filter(n > 5) , aes(x=reason.simplified, y=n)) + geom_bar(stat="identity") + coord_flip() + xlab("Reason for exclusion") + theme(text = element_text(size=20))

## frequency of study
year.rate <- evidence %>% group_by(Publication.Year) %>% summarize(n=length(Publication.Year))

ggplot(tail(year.rate,30)) + geom_bar(aes(x=Publication.Year, y=n), stat="identity") + ylab("number of published studies") +xlab("year published") +theme(text = element_text(size=16))

Prisma report

## total number of papers found
nrow(evidence)
## [1] 3488
## number of papers found outside of WoS
other <- read.csv("data/synthesisdata//other.sources.csv")
nrow(other)
## [1] 0
## number of articles excluded
excludes <- evidence %>% filter(exclude=="yes")
nrow(excludes)
## [1] 3149
## relevant papers
review <- evidence %>% filter(exclude!="yes")
nrow(review)
## [1] 339
## papers for meta
datasets <- read.csv("data//binary.simplified.csv")
meta <- length(unique(datasets$uniqueID))
meta
## [1] 239
prisma(found = 3485,
       found_other = 4,
       no_dupes = 3489,
       screened = 3489,
       screen_exclusions = 0,
       full_text = 3489,
       full_text_exclusions = 3250,
       qualitative = 239, 
       quantitative = 90,
       width = 800, height = 800)

3. Synthesis - Abundance

meta <- read.csv("data//binary.simplified.csv")


## Load packages and functions
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(metafor)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loading 'metafor' package (version 2.1-0). For an overview 
## and introduction to the package please type: help(metafor).
source("meta.evalSE.r") ## Multiple aggregate


## Simplify the estimates column
meta2 <- meta
est <- read.csv("data//Unique.Estimates.Column.csv")
meta2 <- merge(meta2, est, by="Estimate")

## drop Lichen, Microbes and fungi
meta2 <- meta2 %>% filter(Functional.Group != "Producer") %>% filter(Higher.Taxa != "Microbial") ## examine animals only

## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta2 <- merge(meta2, ms.data, by="uniqueID")
## subset for abundance only
meta2 <- meta2 %>% filter(est.reduced == "diversity")
meta2[,"functional.taxa"] <- paste0(meta2$Higher.Taxa, meta2$Functional.Group)

## Create Unique identifier column
meta2[,"UniqueSite"] <- paste(meta2$uniqueID, meta2$Higher.Taxa, meta2$Taxa, meta2$functional.taxa, meta2$estimate.simplified, sep="-")

## convert se to sd
meta2[meta2$Stat=="sd","Value"] <- meta2[meta2$Stat=="sd","Value"] / sqrt(meta2[meta2$Stat=="sd","replicate"] )
meta2[meta2$Stat=="sd","Stat"] <- "se"

## drop comparisons that are not pairwise
meta2 <-  meta2 %>% filter(grazing.reduced == "ungrazed" | grazing.reduced == "grazed")
meta2 <- meta2 %>% filter(!uniqueID %in% c("30","416")) %>% ## drop study 30 and 416 because anomalous 
  filter(uniqueID != "111") ## data is missing in study 111

## Drop anomalous studies
meta2 <- meta2 <- meta2 %>% filter(!uniqueID %in% c("654", "1660")) ## native grazer not ungulate

## Test one higher taxa at a time
meta2 <- meta2 %>% filter(Functional.Group != "Other")

## Test domestic grazers only
meta2 <- meta2 %>%  filter(grazer.simplified %in% c("domestic","both")) %>% filter(!grazer.spp.reduced == "domestic")

## Use function to extract summary statistics for comparisons
## meta.eval  arguments are (meta.data, compare, ids , stats)
grazed.compare <- meta.eval(meta2, grazing.reduced, UniqueSite, Stat)

## Combine the lists into same dataframe
## Rename Columns in second dataframe
grazed.stat <- grazed.compare [[2]] ## extracted statistics 
names(grazed.stat) <- c("UniqueSite","grazed_mean","grazed_se","ungrazed_mean","ungrazed_se","grazed_n","ungrazed_n") ## rename columns to match
grazed.raw <- grazed.compare[[1]] ## calculated statistics from raw values

## Join two dataframes
meta.stat <- rbind(grazed.raw, grazed.stat[, names(grazed.raw)])

## Calculate effect size from SE rather than SD. To calculate variance treat replicates as 1. 
## https://stats.stackexchange.com/questions/152172/meta-analysis-in-r-with-standard-errors-instead-of-standard-deviations-metafor 
meta.stat[,"dummyN"] <- 1
meta.ready <- escalc(n1i = dummyN, n2i = dummyN, m1i = ungrazed_mean, m2i = grazed_mean, sd1i = ungrazed_se, sd2i = grazed_se, data = meta.stat, measure = "ROM", append = TRUE)
  
  ## clean up meta.ready
  meta.ready <- na.omit(meta.ready) ## drop NAs

  
  ## separate out the identifiers
  siteID <- matrix(unlist(strsplit(meta.ready$UniqueSite,"-")),ncol=5, byrow=TRUE)
  siteID <- data.frame(siteID) ## recreate as dataframe
  colnames(siteID) <- c("Study","higher.taxa","taxa","FG","measure") ## add column names
  meta.ready <- cbind(data.frame(meta.ready), siteID)
  
  #random-effects meta-analysis for grazed vs ungrazed plots
  m1 <- rma(yi=yi, vi=vi,  data = meta.ready)
  summary(m1) 
## 
## Random-Effects Model (k = 29; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -4.8711    9.7422   13.7422   16.4066   14.2222   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0427 (SE = 0.0180)
## tau (square root of estimated tau^2 value):      0.2065
## I^2 (total heterogeneity / total variability):   72.39%
## H^2 (total variability / sampling variability):  3.62
## 
## Test for Heterogeneity:
## Q(df = 28) = 109.0878, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.1807  0.0504  3.5848  0.0003  0.0819  0.2794  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  #mixed-effects meta-analysis for grazed vs ungrazed plots
  m2 <- rma(yi=yi, vi=vi, mods=~ -1 + FG,  data = subset(meta.ready, vi<4), digits=4 )
  summary(m2) 
## 
## Mixed-Effects Model (k = 29; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -1.2542    2.5085   22.5085   32.4658   46.9529   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0248 (SE = 0.0151)
## tau (square root of estimated tau^2 value):             0.1574
## I^2 (residual heterogeneity / unaccounted variability): 56.54%
## H^2 (unaccounted variability / sampling variability):   2.30
## 
## Test for Residual Heterogeneity:
## QE(df = 20) = 41.3645, p-val = 0.0033
## 
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 32.6884, p-val = 0.0002
## 
## Model Results:
## 
##                             estimate      se     zval    pval    ci.lb   ci.ub 
## FGInvertebrateDetritivores    0.1089  0.0963   1.1306  0.2582  -0.0799  0.2978 
## FGInvertebrateHerbivorous     0.3712  0.1242   2.9887  0.0028   0.1278  0.6147 
## FGInvertebrateOmnivorous      0.0803  0.2591   0.3098  0.7567  -0.4276  0.5882 
## FGInvertebrateParasitic       0.3580  0.2036   1.7584  0.0787  -0.0410  0.7571 
## FGInvertebratePollinator      0.3032  0.0772   3.9255  <.0001   0.1518  0.4545 
## FGInvertebratePredator        0.2389  0.1384   1.7263  0.0843  -0.0323  0.5101 
## FGVertebrateHerbivorous      -0.0473  0.1315  -0.3597  0.7190  -0.3050  0.2104 
## FGVertebrateOmnivorous       -0.0864  0.1175  -0.7349  0.4624  -0.3167  0.1440 
## FGVertebratePredator         -0.4832  1.0064  -0.4801  0.6312  -2.4557  1.4893 
##  
## FGInvertebrateDetritivores 
## FGInvertebrateHerbivorous    ** 
## FGInvertebrateOmnivorous 
## FGInvertebrateParasitic       . 
## FGInvertebratePollinator    *** 
## FGInvertebratePredator        . 
## FGVertebrateHerbivorous 
## FGVertebrateOmnivorous 
## FGVertebratePredator 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  #The benchmark values for I2 are 25, 50 and 75% for low, moderate and high heterogeneity, respectively (Higgins et al. 2003).
  
  
  ## number of studies per comparison
  nstudies <- meta.ready %>% group_by(FG) %>% summarise(n=length(FG))
  nstudies
## # A tibble: 9 x 2
##   FG                           n
##   <fct>                    <int>
## 1 InvertebrateDetritivores     5
## 2 InvertebrateHerbivorous      4
## 3 InvertebrateOmnivorous       1
## 4 InvertebrateParasitic        1
## 5 InvertebratePollinator       7
## 6 InvertebratePredator         2
## 7 VertebrateHerbivorous        5
## 8 VertebrateOmnivorous         3
## 9 VertebratePredator           1
  #write.csv(meta.ready, "data//effectSize//div.animal.csv", row.names=FALSE)

## compare inclusion of moderators
(.9539-.9352)/.9352 ## explains an additional 2%
## [1] 0.01999572
## Produce a forest plot to determine the effect sizes for each study
forest(m2)

regtest(m2)
## 
## Regression Test for Funnel Plot Asymmetry
## 
## model:     mixed-effects meta-regression model
## predictor: standard error
## 
## test for funnel plot asymmetry: z = 0.7627, p = 0.4456
confint(m1)
## 
##        estimate   ci.lb   ci.ub 
## tau^2    0.0427  0.0149  0.0800 
## tau      0.2065  0.1220  0.2829 
## I^2(%)  72.3950 47.7645 83.1039 
## H^2      3.6225  1.9144  5.9185
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2, ylim=c(1.2,0))

## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)
## 
## Fail-safe N Calculation Using the Rosenthal Approach
## 
## Observed Significance Level: <.0001
## Target Significance Level:   0.05
## 
## Fail-safe N: 384
# ## generate plot with spaces inbetween
# forest(m1, atransf=exp, cex=0.75, ylim=c(-1, 24),
#        order=order(meta.ready$GI.compare,meta.ready$taxa), rows=c(3:4,7,10:16,19:21),
# #         mlab="RE model for all studies", psize=1, slab= paste(meta.ready$Study, meta.ready$taxa, meta.ready$measure))
# 
# addpoly(res.w, row=18, cex=0.75, atransf=exp, mlab="RE model for green wall")
# addpoly(res.r, row= 9, cex=0.75, atransf=exp, mlab="RE model for green roof")
# addpoly(res.rd, row= 6, cex=0.75, atransf=exp, mlab="RE model for roadsides")
# addpoly(res.p, row= 2, cex=0.75, atransf=exp, mlab="RE model for retention ponds")

Create Appendix A

meta <- read.csv("data//binary.simplified.csv")


## Load packages and functions
library(reshape2)
library(metafor)
source("meta.evalSE.r") ## Multiple aggregate


## Simplify the estimates column
meta2 <- meta
est <- read.csv("data//Unique.Estimates.Column.csv")
meta2 <- merge(meta2, est, by="Estimate")

## drop duplicates
ests <- meta2[!duplicated(meta2[,c(1,21:22)]),c(1,21:22)]
ests.used <- ests %>% filter(!est.reduced %in% c("omit", "behavioural","fitness")) ## filter only the variables used

## write estimates
write.csv(ests.used, "data//appendixA.csv", row.names = F)

Animal comparisons

meta.ready <- read.csv("abundanceEff.csv")
  
  #random-effects meta-analysis for grazed vs ungrazed plots
  m1 <- rma(yi=yi, vi=vi,  data = meta.ready, test="t" )
  summary(m1) 
## 
## Random-Effects Model (k = 122; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -140.5963   281.1926   285.1926   290.7842   285.2943   
## 
## tau^2 (estimated amount of total heterogeneity): 0.4202 (SE = 0.0719)
## tau (square root of estimated tau^2 value):      0.6482
## I^2 (total heterogeneity / total variability):   92.78%
## H^2 (total variability / sampling variability):  13.85
## 
## Test for Heterogeneity:
## Q(df = 121) = 2302.4568, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval    pval   ci.lb   ci.ub 
##   0.2089  0.0692  3.0203  0.0031  0.0720  0.3459  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  #mixed-effects meta-analysis for grazed vs ungrazed plots
  m2 <- rma(yi=yi, vi=vi, mods=~ -1 + FG,  data = subset(meta.ready, vi<4), digits=4, test="t" )
  summary(m2) 
## 
## Mixed-Effects Model (k = 122; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -127.2609   254.5219   274.5219   301.7958   276.6788   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.3838 (SE = 0.0696)
## tau (square root of estimated tau^2 value):             0.6195
## I^2 (residual heterogeneity / unaccounted variability): 89.86%
## H^2 (unaccounted variability / sampling variability):   9.86
## 
## Test for Residual Heterogeneity:
## QE(df = 113) = 1723.5800, p-val < .0001
## 
## Test of Moderators (coefficients 1:9):
## F(df1 = 9, df2 = 113) = 3.0013, p-val = 0.0030
## 
## Model Results:
## 
##                             estimate      se     tval    pval    ci.lb    ci.ub 
## FGInvertebrateDetritivores   -0.3829  0.1930  -1.9839  0.0497  -0.7654  -0.0005 
## FGInvertebrateHerbivorous    -0.0423  0.1854  -0.2284  0.8198  -0.4097   0.3250 
## FGInvertebrateOmnivorous      0.3817  0.3808   1.0025  0.3183  -0.3727   1.1362 
## FGInvertebrateParasitic       0.1809  0.2729   0.6630  0.5087  -0.3597   0.7215 
## FGInvertebratePollinator      0.2097  0.2163   0.9697  0.3343  -0.2188   0.6382 
## FGInvertebratePredator        0.1892  0.2211   0.8557  0.3940  -0.2489   0.6273 
## FGVertebrateHerbivorous       0.4626  0.1370   3.3771  0.0010   0.1912   0.7340 
## FGVertebrateOmnivorous        0.2931  0.1947   1.5055  0.1350  -0.0926   0.6788 
## FGVertebratePredator          0.4841  0.1939   2.4971  0.0140   0.1000   0.8682 
##  
## FGInvertebrateDetritivores   * 
## FGInvertebrateHerbivorous 
## FGInvertebrateOmnivorous 
## FGInvertebrateParasitic 
## FGInvertebratePollinator 
## FGInvertebratePredator 
## FGVertebrateHerbivorous     ** 
## FGVertebrateOmnivorous 
## FGVertebratePredator         * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  library(snow)
  hetero <- gosh(m1, subsets=10000, parallel = "snow", ncpus=40)
## Fitting 10000 models (based on random subsets).
  plot(hetero, breaks=100)

### Bias tests  
regtest(m2)
## 
## Regression Test for Funnel Plot Asymmetry
## 
## model:     mixed-effects meta-regression model
## predictor: standard error
## 
## test for funnel plot asymmetry: t = -1.4297, df = 112, p = 0.1556
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2, ylim=c(1.2,0))

## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)
## 
## Fail-safe N Calculation Using the Rosenthal Approach
## 
## Observed Significance Level: <.0001
## Target Significance Level:   0.05
## 
## Fail-safe N: 8156
  meta.ready <- read.csv("diversityEff.csv")
  
  #random-effects meta-analysis for grazed vs ungrazed plots
  m1 <- rma(yi=yi, vi=vi,  data = meta.ready)
  summary(m1) 
## 
## Random-Effects Model (k = 29; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -4.8711    9.7422   13.7422   16.4066   14.2222   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0427 (SE = 0.0180)
## tau (square root of estimated tau^2 value):      0.2065
## I^2 (total heterogeneity / total variability):   72.39%
## H^2 (total variability / sampling variability):  3.62
## 
## Test for Heterogeneity:
## Q(df = 28) = 109.0878, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.1807  0.0504  3.5848  0.0003  0.0819  0.2794  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  hetero <- gosh(m1, subsets=10000, parallel = "snow", ncpus=40)
## Fitting 10000 models (based on random subsets).
  plot(hetero, breaks=100)

  #mixed-effects meta-analysis for grazed vs ungrazed plots
  m2 <- rma(yi=yi, vi=vi, mods=~ -1 + FG,  data = subset(meta.ready, vi<4), digits=4, test="t"  )
  summary(m2) 
## 
## Mixed-Effects Model (k = 29; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -1.2542    2.5085   22.5085   32.4658   46.9529   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0248 (SE = 0.0151)
## tau (square root of estimated tau^2 value):             0.1574
## I^2 (residual heterogeneity / unaccounted variability): 56.54%
## H^2 (unaccounted variability / sampling variability):   2.30
## 
## Test for Residual Heterogeneity:
## QE(df = 20) = 41.3645, p-val = 0.0033
## 
## Test of Moderators (coefficients 1:9):
## F(df1 = 9, df2 = 20) = 3.6320, p-val = 0.0078
## 
## Model Results:
## 
##                             estimate      se     tval    pval    ci.lb   ci.ub 
## FGInvertebrateDetritivores    0.1089  0.0963   1.1306  0.2716  -0.0920  0.3099 
## FGInvertebrateHerbivorous     0.3712  0.1242   2.9887  0.0073   0.1121  0.6303 
## FGInvertebrateOmnivorous      0.0803  0.2591   0.3098  0.7599  -0.4603  0.6208 
## FGInvertebrateParasitic       0.3580  0.2036   1.7584  0.0940  -0.0667  0.7828 
## FGInvertebratePollinator      0.3032  0.0772   3.9255  0.0008   0.1421  0.4643 
## FGInvertebratePredator        0.2389  0.1384   1.7263  0.0997  -0.0498  0.5276 
## FGVertebrateHerbivorous      -0.0473  0.1315  -0.3597  0.7228  -0.3216  0.2270 
## FGVertebrateOmnivorous       -0.0864  0.1175  -0.7349  0.4709  -0.3316  0.1588 
## FGVertebratePredator         -0.4832  1.0064  -0.4801  0.6364  -2.5825  1.6161 
##  
## FGInvertebrateDetritivores 
## FGInvertebrateHerbivorous    ** 
## FGInvertebrateOmnivorous 
## FGInvertebrateParasitic       . 
## FGInvertebratePollinator    *** 
## FGInvertebratePredator        . 
## FGVertebrateHerbivorous 
## FGVertebrateOmnivorous 
## FGVertebratePredator 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ### Bias tests  
regtest(m2)
## 
## Regression Test for Funnel Plot Asymmetry
## 
## model:     mixed-effects meta-regression model
## predictor: standard error
## 
## test for funnel plot asymmetry: t = 0.7627, df = 19, p = 0.4550
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2, ylim=c(0.8,0))

## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)
## 
## Fail-safe N Calculation Using the Rosenthal Approach
## 
## Observed Significance Level: <.0001
## Target Significance Level:   0.05
## 
## Fail-safe N: 384

Testing the effects of grazing on producers

meta <- read.csv("data//plantAnalysis.csv")
meta2 <- meta
  
## Load packages and functions
library(reshape2)
library(metafor)
source("meta.evalSE.r") ## Multiple aggregate
  
  
  
  ## join manuscript data
  ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
  meta2 <- merge(meta2, ms.data, by="uniqueID")

  ## Drop mychoriza because only two studies
  
    ## Create Unique identifier column
  meta2 <- meta2 %>% filter(estimate.simplified != "omit")
  meta2[,"UniqueSite"] <- paste(meta2$uniqueID, meta2$Taxa, meta2$estimate.simplified, sep="-")
  
  ## convert se to sd
   meta2[meta2$Stat=="sd","Value"] <- meta2[meta2$Stat=="sd","Value"] / sqrt(meta2[meta2$Stat=="sd","replicate"] )
   meta2[meta2$Stat=="sd","Stat"] <- "se"
  
  ## drop comparisons that are not pairwise
  meta2 <-  meta2 %>% filter(grazing.reduced == "ungrazed" | grazing.reduced == "grazed")
  meta2 <- meta2 %>% filter(uniqueID != "30") ## drop study 30 because anomalous 
  
  ## Drop anomalous studies
  meta2 <- meta2 <- meta2 %>% filter(!uniqueID %in% c("654", "1660")) ## native grazer not ungulate
  meta2 <- meta2 %>% filter(estimate.simplified != "fitness") ## drop fitness because only study 

  ## Test domestic grazers only
  meta2 <- meta2 %>% filter(grazer.simplified %in% c("domestic","both")) %>% 
    filter(grazer.spp.reduced != "domestic") ## select only sheep and cattle comparisons

    ## Use function to extract summary statistics for comparisons
  ## meta.eval  arguments are (meta.data, compare, ids , stats)
  grazed.compare <- meta.eval(meta2, grazing.reduced, UniqueSite, Stat)
  
  ## Combine the lists into same dataframe
  ## Rename Columns in second dataframe
  grazed.stat <- grazed.compare [[2]] ## extracted statistics 
  names(grazed.stat) <- c("UniqueSite","grazed_mean","grazed_se","ungrazed_mean","ungrazed_se","grazed_n","ungrazed_n") ## rename columns to match
  grazed.raw <- grazed.compare[[1]] ## calculated statistics from raw values
  
  ## Join two dataframes
  meta.stat <- rbind(grazed.raw, grazed.stat[, names(grazed.raw)])
  
  
## Calculate effect size from SE rather than SD. To calculate variance treat replicates as 1. 
## https://stats.stackexchange.com/questions/152172/meta-analysis-in-r-with-standard-errors-instead-of-standard-deviations-metafor 
meta.stat[,"dummyN"] <- 1
meta.ready <- escalc(n1i = dummyN, n2i = dummyN, m1i = ungrazed_mean, m2i = grazed_mean, sd1i = ungrazed_se, sd2i = grazed_se, data = meta.stat, measure = "ROM", append = TRUE)
## Warning in escalc(n1i = dummyN, n2i = dummyN, m1i = ungrazed_mean, m2i =
## grazed_mean, : Some 'yi' and/or 'vi' values equal to +-Inf. Recoded to NAs.
  ## clean up meta.ready
  meta.ready <- na.omit(meta.ready) ## drop NAs

  
  ## separate out the identifiers
  siteID <- matrix(unlist(strsplit(meta.ready$UniqueSite,"-")),ncol=3, byrow=TRUE)
  siteID <- data.frame(siteID) ## recreate as dataframe
  colnames(siteID) <- c("Study","taxa","measure") ## add column names
  meta.ready <- cbind(data.frame(meta.ready), siteID)
  
  #random-effects meta-analysis for grazed vs ungrazed plots
  m1 <- rma(yi=yi, vi=vi,  data = meta.ready)
## Warning in rma(yi = yi, vi = vi, data = meta.ready): There are outcomes with
## non-positive sampling variances.
## Warning in rma(yi = yi, vi = vi, data = meta.ready): Cannot compute Q-test, I^2,
## or H^2 with non-positive sampling variances.
  summary(m1) 
## 
## Random-Effects Model (k = 65; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -60.9895  121.9791  125.9791  130.2968  126.1758   
## 
## tau^2 (estimated amount of total heterogeneity): 0.3054 (SE = 0.0621)
## tau (square root of estimated tau^2 value):      0.5527
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.2942  0.0741  3.9692  <.0001  0.1489  0.4395  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ## get last year grazed as covariate
  yr.grazed <- meta2 %>% group_by(uniqueID) %>% summarize(yr=max(yr.grazed))
  colnames(yr.grazed)[1] <- "Study"
  site.yr <- merge(siteID, yr.grazed, by="Study")
  meta.ready[,"yr.grazed"] <- site.yr$yr
  meta.ready[,"previously.grazed"] <- as.factor(ifelse(is.na(meta.ready$yr.grazed),0,1))
  
  #mixed-effects meta-analysis for grazed vs ungrazed plots
  m2 <- rma(yi=yi, vi=vi, mods=~ -1 + measure ,  data = subset(meta.ready, vi!=0))
  summary(m2) 
## 
## Mixed-Effects Model (k = 61; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -50.3340  100.6679  110.6679  120.8832  111.8444   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.2471 (SE = 0.0553)
## tau (square root of estimated tau^2 value):             0.4970
## I^2 (residual heterogeneity / unaccounted variability): 95.65%
## H^2 (unaccounted variability / sampling variability):   23.00
## 
## Test for Residual Heterogeneity:
## QE(df = 57) = 485.3761, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 36.9574, p-val < .0001
## 
## Model Results:
## 
##                            estimate      se     zval    pval    ci.lb   ci.ub 
## measureabundance             0.4267  0.0872   4.8932  <.0001   0.2558  0.5976 
## measuredebris                0.9367  0.2710   3.4568  0.0005   0.4056  1.4678 
## measurediversity            -0.1251  0.1415  -0.8839  0.3768  -0.4024  0.1523 
## measureMycorrhizaColonize    0.1887  0.3545   0.5324  0.5944  -0.5060  0.8835 
##  
## measureabundance           *** 
## measuredebris              *** 
## measurediversity 
## measureMycorrhizaColonize 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  ## number of studies
  meta.ready %>% group_by(measure) %>% summarize(n=length(measure))
## # A tibble: 4 x 2
##   measure                n
##   <fct>              <int>
## 1 abundance             41
## 2 debris                 5
## 3 diversity             17
## 4 MycorrhizaColonize     2
  fsn(yi, vi, data=meta.ready)
## 
## Fail-safe N Calculation Using the Rosenthal Approach
## 
## Observed Significance Level: NA
## Target Significance Level:   0.05
## 
## Fail-safe N: NaN
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2)

regtest(m2)
## 
## Regression Test for Funnel Plot Asymmetry
## 
## model:     mixed-effects meta-regression model
## predictor: standard error
## 
## test for funnel plot asymmetry: z = 0.9919, p = 0.3212
## write
write.csv(subset(meta.ready, vi<4), "plantEff.csv", row.names=F)

Plant analysis and GOSH plots

  meta.ready <- read.csv("plantEff.csv")
  
  ## Abundance
  #random-effects meta-analysis for grazed vs ungrazed plots
  m1 <- rma(yi=yi, vi=vi,  data = subset(meta.ready, measure == "abundance" & vi!=0))
  summary(m1) 
## 
## Random-Effects Model (k = 40; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -32.6182   65.2364   69.2364   72.5635   69.5697   
## 
## tau^2 (estimated amount of total heterogeneity): 0.2399 (SE = 0.0659)
## tau (square root of estimated tau^2 value):      0.4898
## I^2 (total heterogeneity / total variability):   96.06%
## H^2 (total variability / sampling variability):  25.38
## 
## Test for Heterogeneity:
## Q(df = 39) = 340.6311, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.4261  0.0861  4.9478  <.0001  0.2573  0.5949  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  hetero <- gosh(m1, subsets=10000, parallel = "snow", ncpus=10)
## Fitting 10000 models (based on random subsets).
  plot(hetero, breaks=100)

  ## Diversity
  #random-effects meta-analysis for grazed vs ungrazed plots
  m2 <- rma(yi=yi, vi=vi,  data = subset(meta.ready, measure == "diversity" & vi!=0))
  summary(m2) 
## 
## Random-Effects Model (k = 14; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -0.9716    1.9432    5.9432    7.0731    7.1432   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0321 (SE = 0.0208)
## tau (square root of estimated tau^2 value):      0.1793
## I^2 (total heterogeneity / total variability):   67.07%
## H^2 (total variability / sampling variability):  3.04
## 
## Test for Heterogeneity:
## Q(df = 13) = 34.4838, p-val = 0.0010
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb   ci.ub 
##  -0.1123  0.0630  -1.7805  0.0750  -0.2358  0.0113  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  hetero <- gosh(m2, subsets=10000, parallel = "snow", ncpus=10)
## Fitting 10000 models (based on random subsets).
  plot(hetero, breaks=100)

Separated Plant responses

Compare regional patterns

## Load in effect sizes
occurdat<-list.files( pattern="Eff.csv$",full=T)
animal.abd <- read.csv(occurdat[2])
# animal.abd <- subset(animal.abd, measure=="diversity" & vi!=0) ## one measure at a time for plants
colnames(animal.abd)[colnames(animal.abd)=="Study"] <- "uniqueID"
 
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta3 <- merge(animal.abd, ms.data, by="uniqueID")
meta3[,c("yi")] <- meta3[,c("yi")] ## inverse effect

library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
# 
# ## Convert lat lon to spatial points
# meta3 <- meta3 %>% filter(lat != "") ## remove blanks for lat and lon
# meta3 <- meta3[!grepl(";", meta3$lat),] ## remove multiple entries
# meta3$lat <- as.numeric(levels(meta3$lat))[meta3$lat] ## convert from factor to number
# meta3$lon <- as.numeric(levels(meta3$lon))[meta3$lon] ## convert from factor to number
# 
# ## generate dataframe
# gps <- data.frame(uniqueID = meta3$uniqueID, lon=meta3$lon, lat=meta3$lat)
# gps <- gps[!duplicated(gps$uniqueID),]
# 
# 
# ## assign spatial points
# crs.world <-CRS("+proj=longlat +datum=WGS84")
# coordinates(gps) <- ~lon+lat
# proj4string(gps) <- crs.world
# 
# ## load rasters
# temp <- raster("C:\\Users\\Fitz\\Downloads\\bioclim2\\bio01.tif")
# prec <- raster("C:\\Users\\Fitz\\Downloads\\bioclim2\\bio12.tif")
# 
# temp2 <- temp
# 
# clim <- stack(temp,prec)
# 
# 
# ## extract temperature and preciptiation
# gps.clim <- data.frame(raster::extract(clim, gps))
# gps.clim["uniqueID"] <- gps$uniqueID


## get CRU data
CRU <- read.csv("data//climateCRU.csv")
CRU <- na.omit(CRU)

## compare climate & grazing duration variables against effect size


## get last year grazed as covariate
meta.ready <- merge(meta3, CRU, by="uniqueID")

### Temperature only
#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~ poly(temp,2),  data = meta.ready , test="t" )
summary(m3) 
## 
## Mixed-Effects Model (k = 27; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -4.1041    8.2082   16.2082   20.9204   18.3135   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0377 (SE = 0.0182)
## tau (square root of estimated tau^2 value):             0.1942
## I^2 (residual heterogeneity / unaccounted variability): 68.63%
## H^2 (unaccounted variability / sampling variability):   3.19
## R^2 (amount of heterogeneity accounted for):            16.62%
## 
## Test for Residual Heterogeneity:
## QE(df = 24) = 68.9193, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 24) = 2.2044, p-val = 0.1322
## 
## Model Results:
## 
##                 estimate      se     tval    pval    ci.lb   ci.ub 
## intrcpt           0.1789  0.0503   3.5571  0.0016   0.0751  0.2827  ** 
## poly(temp, 2)1   -0.2404  0.2547  -0.9435  0.3548  -0.7661  0.2854     
## poly(temp, 2)2   -0.4595  0.2378  -1.9320  0.0653  -0.9504  0.0314   . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## line
preds <- predict(m3, temp=meta.ready$temp)

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")

## generate column for FG
meta.ready[,"trophic"] <- gsub("^.*?ertebrate","",meta.ready$FG)

### Plot against temperature
ggplot(meta.ready, aes(x=temp, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa), size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Mean annual temperature (°C)") + ylab("Effect of livestock exclusion (LRR)") + geom_line(data=data.frame(preds), aes(x=meta.ready$temp, y=pred), lwd=2) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$temp, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + scale_color_manual(values=cbPalette)+xlim(0,20) + annotate("text", x=2, y=-.6, label="animal \n diversity", size=6)
## Warning: Ignoring unknown aesthetics: y

### Preciptiation

test <- subset(meta.ready, FG == "VertebrateHerbivorous" | FG == "InvertebrateHerbivorous")

#mixed-effects meta-analysis for grazed vs ungrazed plots
m4 <- rma(yi=yi, vi=vi, mods=~ poly(precip,2),  data = subset(meta.ready))
summary(m4) 
## 
## Mixed-Effects Model (k = 27; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -5.2128   10.4255   18.4255   23.1377   20.5308   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0448 (SE = 0.0203)
## tau (square root of estimated tau^2 value):             0.2116
## I^2 (residual heterogeneity / unaccounted variability): 73.22%
## H^2 (unaccounted variability / sampling variability):   3.73
## R^2 (amount of heterogeneity accounted for):            1.05%
## 
## Test for Residual Heterogeneity:
## QE(df = 24) = 86.7515, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 1.4233, p-val = 0.4908
## 
## Model Results:
## 
##                   estimate      se     zval    pval    ci.lb   ci.ub 
## intrcpt             0.1685  0.0537   3.1358  0.0017   0.0632  0.2738  ** 
## poly(precip, 2)1    0.1461  0.2857   0.5112  0.6092  -0.4139  0.7060     
## poly(precip, 2)2   -0.2959  0.2836  -1.0433  0.2968  -0.8517  0.2600     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## line
preds <- predict(m4, precip=meta.ready$precip)


### Plot against temperature
ggplot(meta.ready, aes(x=precip, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa),size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Mean annual preciptiation (mm)") + ylab("Effect of livestock exclusion (LRR)") + scale_color_manual(values=cbPalette) + 
  geom_line(data=data.frame(preds), aes(x=meta.ready$precip, y=pred), lwd=2) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$precip, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) 
## Warning: Ignoring unknown aesthetics: y

## Time since grazed Grazed

## Load in effect sizes
occurdat<-list.files( pattern="Eff.csv$",full=T)
animal.abd <- read.csv(occurdat[1])
colnames(animal.abd)[colnames(animal.abd)=="Study"] <- "uniqueID"
 
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta3 <- merge(animal.abd, ms.data, by="uniqueID")
meta3$last.grazing.event <- as.numeric(meta3$last.grazing.event)
meta3[,c("yi")] <- meta3[,c("yi")]
meta3[,"yr.grazed"] <- meta3$last.grazing.event

## compare climate & grazing duration variables against effect size

## get last year grazed as covariate
meta.ready <- subset(meta3, yr.grazed>0& yr.grazed < 60)


#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~  poly(yr.grazed,1),  data = meta.ready, test="t")
summary(m3) 
## 
## Mixed-Effects Model (k = 107; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -113.7611   227.5222   233.5222   241.4841   233.7598   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.2984 (SE = 0.0593)
## tau (square root of estimated tau^2 value):             0.5462
## I^2 (residual heterogeneity / unaccounted variability): 86.80%
## H^2 (unaccounted variability / sampling variability):   7.58
## R^2 (amount of heterogeneity accounted for):            9.95%
## 
## Test for Residual Heterogeneity:
## QE(df = 105) = 620.3301, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 105) = 7.0875, p-val = 0.0090
## 
## Model Results:
## 
##                     estimate      se    tval    pval   ci.lb   ci.ub 
## intrcpt               0.1458  0.0652  2.2352  0.0275  0.0165  0.2751   * 
## poly(yr.grazed, 1)    1.8068  0.6787  2.6622  0.0090  0.4611  3.1524  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## generate column for FG
meta.ready[,"trophic"] <- gsub("^.*?ertebrate","",meta.ready$FG)


## line
preds <- predict(m3, yr.grazed=meta.ready$yr.grazed)

# cbPalette <- c("#77564C", "#4381C1", "#9AD4D6", "#4FB286", "#564787", "#08415C")

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")


## animals
plot1 <- ggplot(meta.ready, aes(x=yr.grazed, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa), size=4, width=1) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Time since exclusion (yrs)") + ylab("Effect of livestock exclusion (LRR)") + geom_line(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred), lwd=1) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + scale_color_manual(values=cbPalette)+ 
  annotate("text", x=32, y=2.4, label="animal abundance", size=8)
## Warning: Ignoring unknown aesthetics: y
#### diversity
## Load in effect sizes
occurdat<-list.files( pattern="Eff.csv$",full=T)
animal.abd <- read.csv(occurdat[2])
colnames(animal.abd)[colnames(animal.abd)=="Study"] <- "uniqueID"
 
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta3 <- merge(animal.abd, ms.data, by="uniqueID")
meta3$last.grazing.event <- as.numeric(meta3$last.grazing.event)
meta3[,c("yi")] <- meta3[,c("yi")]
meta3[,"yr.grazed"] <- meta3$last.grazing.event

## compare climate & grazing duration variables against effect size

## get last year grazed as covariate
meta.ready <- subset(meta3, yr.grazed>0& yr.grazed < 60)


#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~  poly(yr.grazed,2),  data = meta.ready, test="t")
summary(m3) 
## 
## Mixed-Effects Model (k = 26; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -1.7467    3.4934   11.4934   16.0353   13.7156   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0340 (SE = 0.0172)
## tau (square root of estimated tau^2 value):             0.1843
## I^2 (residual heterogeneity / unaccounted variability): 64.97%
## H^2 (unaccounted variability / sampling variability):   2.85
## R^2 (amount of heterogeneity accounted for):            13.85%
## 
## Test for Residual Heterogeneity:
## QE(df = 23) = 64.9276, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 23) = 2.0041, p-val = 0.1576
## 
## Model Results:
## 
##                      estimate      se     tval    pval    ci.lb   ci.ub 
## intrcpt                0.1981  0.0492   4.0299  0.0005   0.0964  0.2999  *** 
## poly(yr.grazed, 2)1   -0.0935  0.2567  -0.3642  0.7191  -0.6244  0.4375      
## poly(yr.grazed, 2)2   -0.5362  0.2679  -2.0020  0.0572  -1.0903  0.0178    . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## generate column for FG
meta.ready[,"trophic"] <- gsub("^.*?ertebrate","",meta.ready$FG)


## line
preds <- predict(m3, yr.grazed=meta.ready$yr.grazed)

# cbPalette <- c("#77564C", "#4381C1", "#9AD4D6", "#4FB286", "#564787", "#08415C")

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")


## animals
plot2 <- ggplot(meta.ready, aes(x=yr.grazed, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa), size=4, width=1) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Time since exclusion (yrs)") + ylab("Effect of livestock exclusion (LRR)") + geom_line(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred), lwd=1) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + scale_color_manual(values=cbPalette)+ 
  annotate("text", x=32, y=1, label="animal diversity", size=8)
## Warning: Ignoring unknown aesthetics: y
### Old figure
# ### calculate predicted risk ratios for full grazing time gradient
# preds <- predict(m3, newmods=c(0:35), transf=exp)
# 
# ### radius of points will be proportional to the inverse standard errors
# ### hence the area of the points will be proportional to inverse variances
# size <- 1 / sqrt(meta.ready$vi)
# size <- size / max(size) *0.95
# 
# ### set up plot (risk ratios on y-axis, absolute temperature on x-axis)
# plot(NA, NA, xlim=c(0,35), ylim=c(0.5,5),
#      xlab="Year since control was last grazed", ylab="Effect of grazing relative to control (LRR)",
#      las=1, bty="l", log="y", cex.lab=1.5, cex.axis=1.2)
# 
# ### add points
# col <- data.frame(FG=unique(meta.ready$FG), colors= c("#999999","#E69F00", "#56B4E9", "#D55E00","#0072B2","#E69F00","#0072B2"))
# ptshape <- data.frame(FG=unique(meta.ready$FG), shape= c(21,21,21,21,21,23,23))
# meta.ready <- merge(meta.ready, col, by="FG")
# meta.ready <- merge(meta.ready, ptshape, by="FG")
# points(meta.ready$yr.grazed, exp(meta.ready$yi), pch=meta.ready$shape, bg=as.character(meta.ready$colors), cex=2)
# arrows(meta.ready$yr.grazed, exp(meta.ready$yi)-meta.ready$vi, ## lower
#        meta.ready$yr.grazed, exp(meta.ready$yi)+meta.ready$vi, ## upper
#        length=0.05, angle=90, code=3)
# 
# # symbols(meta.ready$yr.grazed, exp(meta.ready$yi), circles=size, inches=FALSE, add=TRUE, bg=as.character(meta.ready$colors))
# 
# legend(0,5.5, legend=col$FG, pt.bg=as.character(col$colors), pch=meta.ready$shape, cex=0.8)
#                    
# ### add predicted values (and corresponding CI bounds)
# lines(0:35, preds$pred, lwd=2)
# lines(0:35, preds$ci.lb, lty="dashed")
# lines(0:35, preds$ci.ub, lty="dashed")

Taxa that are tested

### Load main data
meta <- read.csv("data//binary.simplified.csv")

## load studies used
div <- read.csv("diversityEff.csv")
abd <- read.csv("abundanceEff.csv")

studiesUsed <- c(div$Study, abd$Study)
dataused <- meta %>% filter(uniqueID %in% studiesUsed)

## count taxa 
taxa <- dataused %>% group_by(Taxa) %>% summarize(n=length(unique(uniqueID)))
taxa
## # A tibble: 17 x 2
##    Taxa             n
##    <fct>        <int>
##  1 ""               4
##  2 "AMF"            1
##  3 "Amphibian"      2
##  4 "Amphibians"     2
##  5 "Annelida"       1
##  6 "Arthropoda"    41
##  7 "Aves"          12
##  8 "Insecta"        1
##  9 "Mammalia"      31
## 10 "Microbes"       1
## 11 "Mollusca"       1
## 12 "Nematoda"       2
## 13 "Nematode"       1
## 14 "Orthoptera"     1
## 15 "plants"         1
## 16 "Plants"        27
## 17 "Reptilia"       8
## unique studies
length(unique(studiesUsed))
## [1] 90
## list species
dataused %>% group_by(Detailed.Taxa ) %>% summarize(n=length(unique(Detailed.Taxa ))) %>%  data.frame(.)
##     Detailed.Taxa n
## 1                 1
## 2           Acari 1
## 3      Amphibians 1
## 4        Annelida 1
## 5       Arachnida 1
## 6         Araneae 1
## 7      Arthropoda 1
## 8             Bee 1
## 9         Bovidae 1
## 10     Coleoptera 1
## 11     Collembola 1
## 12 Dasyuromorphia 1
## 13        Diptera 1
## 14        Equidae 1
## 15     Formicidae 1
## 16     Gastropoda 1
## 17      Hemiptera 1
## 18    Hymenoptera 1
## 19       Isoptera 1
## 20     Lagomorpha 1
## 21    Lepidoptera 1
## 22     Mycorrhiza 1
## 23     Orthoptera 1
## 24          plant 1
## 25     Psocoptera 1
## 26        Rodenta 1
## 27       Squamata 1
## 28   Thysanoptera 1

Cascade effects

## load studies used
div <- read.csv("diversityEff.csv")
abd <- read.csv("abundanceEff.csv")
plt <- read.csv("plantEff.csv")

## merge data
plt <- subset(plt, measure=="diversity")

## select matching studies
abd <- abd[abd$Study %in% plt$Study,]
plt <- plt[plt$Study %in% abd$Study,]

## generate dataframe
plt <- data.frame(Study=plt$Study, plteff = plt$yi)
pltabd <- merge(plt, abd, by="Study", all.x=T)

## Model
m1 <- rma(yi=yi, vi=vi, mods=~ plteff,  data = pltabd, test="t")
summary(m1)
## 
## Mixed-Effects Model (k = 15; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##  -6.5583   13.1166   19.1166   20.8114   21.7832   
## 
## tau^2 (estimated amount of residual heterogeneity):     0 (SE = 0.0234)
## tau (square root of estimated tau^2 value):             0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability):   1.00
## R^2 (amount of heterogeneity accounted for):            100.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 13) = 12.5828, p-val = 0.4805
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 13) = 13.2080, p-val = 0.0030
## 
## Model Results:
## 
##          estimate      se    tval    pval   ci.lb   ci.ub 
## intrcpt    0.3604  0.0730  4.9378  0.0003  0.2027  0.5181  *** 
## plteff     1.1351  0.3123  3.6343  0.0030  0.4603  1.8098   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## separate trophic group
pltabd[,"trophic"] <- gsub("^.*?ertebrate","",pltabd$FG)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")

## predict line
preds <- predict(m1, newdata=data.frame(plteff=pltabd$plteff))

##plot
ggplot(data=pltabd, aes(x=plteff, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+geom_vline(xintercept=0, lwd=1, lty=2)+ geom_point(aes(color=trophic, shape=higher.taxa), size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Effect of livestock exclusion on plant diversity (LRR)") + ylab("Effect of livestock exclusion on animal abundance (LRR)") + geom_line(data=data.frame(preds), aes(x=pltabd$plteff, y=pred), lwd=1) +
geom_ribbon(data=data.frame(preds), aes(x=pltabd$plteff, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) +
  scale_color_manual(values=cbPalette) 
## Warning: Ignoring unknown aesthetics: y

  # annotate("text", x=-.25, y=-1.6, label="plant diversity \n animal abundance", size=6)


### Animal diversity
## merge data

## select matching studies
div <- div[div$Study %in% plt$Study,]
plt <- plt[plt$Study %in% div$Study,]

## generate dataframe
#plt <- data.frame(Study=plt$Study, plteff = plt$yi)
pltdiv <- merge(plt, div, by="Study", all.x=T)

## Model
m1 <- rma(yi=yi, vi=vi, mods=~  poly(plteff,1),  data = pltdiv, test="t")
summary(m1)
## 
## Mixed-Effects Model (k = 8; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##   1.9345   -3.8691    2.1309    1.5062   14.1309   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0098 (SE = 0.0158)
## tau (square root of estimated tau^2 value):             0.0992
## I^2 (residual heterogeneity / unaccounted variability): 35.87%
## H^2 (unaccounted variability / sampling variability):   1.56
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 6) = 8.6282, p-val = 0.1956
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 6) = 0.0059, p-val = 0.9413
## 
## Model Results:
## 
##                  estimate      se     tval    pval    ci.lb   ci.ub 
## intrcpt            0.2353  0.0604   3.8940  0.0080   0.0874  0.3831  ** 
## poly(plteff, 1)   -0.0117  0.1524  -0.0768  0.9413  -0.3847  0.3612     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## separate trophic group
pltdiv[,"trophic"] <- gsub("^.*?ertebrate","",pltdiv$FG)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")

## predict line
preds <- predict(m1, plteff=pltdiv$plteff)

##plot
ggplot(data=pltdiv, aes(x=plteff, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2) + geom_vline(xintercept=0, lwd=1, lty=2)+  geom_point(aes(color=trophic, shape=higher.taxa), size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Effect of livestock on plants (LRR)") + ylab("Effect of livestock on animals (LRR)") +# geom_line(data=data.frame(preds), aes(x=pltdiv$plteff, y=pred), lwd=1) +
#geom_ribbon(data=data.frame(preds), aes(x=pltdiv$plteff, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + 
  scale_color_manual(values=cbPalette)+ ylim(-0.5,0.25)+
  annotate("text", x=-.3, y=-.4, label="plant diversity \n animal diversity", size=6)
## Warning: Removed 3 rows containing missing values (geom_point).